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Abstract. We present a system of interpolating splines with first-order and 
approximate second-order geometric continuity. The curves are easily com- 
puted in linear time by solving a diagonally dominant, tridiagonal system of 
linear equations. Emphasis is placed on the need to find aesthetically pleasing 
curves in a wide range of circumstances; favorable results are obtained even 
when the knots are very unequally spaced or widely separated. The curves are 
invariant under translation, rotation, and scaling, and the effects of a local 
change fall off exponentially as one moves away from the disturbed knot. 

Approximate second-order continuity is achieved by using a linear “mock 
curvature” function in place of the actual endpoint curvature for each spline 
segment and choosing tangent directions at knots so as to equalize these. This 
avoids extraneous solutions and other forms of undesirable behavior without 
seriously compromising the quality of the results. 

The actual spline segments can come from any family of curves whose 
endpoint curvatures can be suitably approximated, but we propose a specific 
family of parametric cubics. There is freedom to allow tangent directions and 
“tension” parameters to be specified at knots, and special “curl” parameters 
may be given for additional control near the endpoints of open curves. 


1. Introduction 


The problem of fitting a smooth curve through a set of points on the plane has 
many important applications in computer graphics, computer-aided design, and 
typesetting. Often there is no preexisting curve to approximate except possibly a 
freehand drawing, and the only requirement is to find an aesthetically pleasing 
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curve that the computer can easily manipulate. For some interactive applications 
the curves can be controlled by manipulating points that do not lie on the curve, 
but many applications require the control points to lie on the curve. For example, 
the control points may be obtained by digitizing key points on a drawing, or there 
may be a priori knowledge that the curve must pass through certain points. 

Suppose the curve must pass through points Zo, Z),...,Z,; and either z) and 
z, are to be the endpoints of the curve or Zp) = z,, and the curve is to be a closed 
loop. Optionally, there may be direction vectors w, specifying the slope of the 
curve at some z,, For example, some of the z, may have been selected as vertical 
extrema so that the curve must pass through them horizontally. It is desirable for 
the curve to be invariant under translation, rotation, and scaling in the sense that 
if T is such a transformation, then applying T to the computed curve should yield 
the same result as computing a new curve through Tz, for 0 <i <n with direction 
vectors Tw, — T(0,0). 

The curve should have at least approximate continuity of slope and curvature 
at knots where no directions are given, and it would also be desirable to have 
some notion of extensibility and locality. A system of splines is extensible if the 
curve generated from knots Zo, z;,...,Z, is identical to that generated from knots 
Z6, Z{,.-., 2441 Where z/ = z, for i < k, z{/=2,_, for i >k, and z; is on the curve 
segment joining z,_, and z,. In other words, adding a new knot already on the 
curve must not change it. In practice it is extremely difficult to achieve exact 
extensibility. The only well-known extensible spline family is the “curve of least 
energy” that minimizes the integral of squared curvature with respect to arc 
length [5, 9], but this curve is difficult to work with. It is interesting to note that 
when the knots are nearly collinear, the curve of least energy approaches the 
simple nonparametric cubic spline passing through the given knots with continu- 
ous second derivative. The splines that we deal with here will share this property. 

The concept of locality is that if one of the knots or direction vectors is 
perturbed, the changes should be confined to a few surrounding spline segments. 
Specifically, there should be some constant k such that changes to z; or w; effect 
only the curve segments between z,_, and z,,,. As the example of Fig. 1 shows, 
it is difficult to have both locality and continuity of curvature even when the 
knots are collinear. If the knots zo, z,, Z2, and z, are as shown in the figure but 
Wy is in the direction of z, — zo, then the desired curve is obviously a straight line. 
if wo is changed, then locality demands that the spline segments between z, and 
n must remain straight, yet there is no way a cubic curve can join such a straight 
line with continuous curvature. 

Most local interpolating splines of moderate complexity do not have continu- 
ous curvature. Perhaps the best known splines of this type are the cardinal splines 
given by a cubic function z(t), where the velocity vector z’(i) at knot z, is 
4(z,;.1— Z;-1)- Figure 2 shows an example of cardinal splines in comparison to 
the splines that will be developed in this paper. The price paid for locality is that 
Fig. 2(b) has wild discontinuities in curvature while Fig. 2(a) does not. 


Pee O Fig. 1. The effect of changing w while preserving 
Zo EA Ze —z; exact locality. 
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(a) (b) (c) 


Fig. 2. (a) Cubics with mock curvature constraints. (b) Cardinal splines. (c) Barsky-DeRose splines 
with B, = ilz, +1 a z\l/iiz, — 2,- ill. 


Other well-known local cubic interpolating splines can be viewed as generali- 
zations of cardinal splines. To see this, it is convenient to parameterize the cubic 
spline in terms of the incoming and outgoing velocity vectors D7 and D;* at each 
knot so that the spline segment between z; and z,,, has Bézier control points z,, 
2Z;41~ Dp’ and z,;,,. For cardinal splines as described above, 
D7 = D} = 3(2;41— 2;-1); in the most general form 


D; = a,(z,-2;-,) + 6(z,41-2;) and 


D* = c(z;— z;-1) + d(Zj44- Zi), (1) 


where a,, b,, c, and d, are parameters that may be used to control the shape of 
the curve. For instance, Kochanek and Bartels [7] give a scheme where the user 
can give continuity, tension, and bias parameters for each knot, and these 
determine a,, b, c,» and d;. The default settings of these parameters result in 
a,=b,=c, = d,;= 4, giving simple cardinal splines as shown in Fig. 2(b). 
Another approach to local interpolating splines can be found in [3], where 
DeRose and Barsky give an infinite family of spline curves based on two integer 
parameters: One parameter gives the desired order of continuity, and the other 
determines whether the splines interpolate the knots or just pass near them. 
DeRose and Barsky give explicit equations for cubic interpolating splines of this 
form with first-order continuity, but their approach requires polynomials of 
degree at least 5 in order to obtain interpolating splines with continuous curva- 
ture. In the cubic interpolating splines derived explicitly in [3], each knot z; has a 
Single shape parameter £, associated with it. The splines are exactly those 
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(a) (b) (c) 


Fig. 3. (a) Natural cubics. (b) Cubics with chordal parameterization. (c) Cubics with mock curvature 
constraints. 


produced by (1) where 


E: SETE A -B j ea 
w= Bt+t 7% BRI FT BAT MS AT BFT 


The default suggested in [3] is to set all 8, =1 so as to obtain cardinal splines as in 
Fig. 2(b), but somewhat better results can be obtained by setting B; equal to the 
ratio of Euclidean lengths |{z,,, ~ z,||/||z; — 2; || as shown in Fig. 2(c). 

In this article we shall see how to obtain smother curves by sacrificing locality 
and settling for an exponential decline in influence of local changes. The best 
known interpolating splines that behave this way are the C? continuous paramet- 
ric cubic splines, sometimes called “natural cubic splines.” If no directions are 
given, there is a unique piecewise parametric cubic, closed curve that is C? 
continuous with respect to the parameter and passes through n given points in 
order. This curve can easily be computed by solving linear equations. 

As Epstein shows in [4], natural cubic splines do not perform well for 
unequally spaced knots because the spacing of parameter values at knots does not 
reflect the spacing of the knots. Better results can be obtained by setting the 
parameter at each knot z; to a value t;, where t; — ¢;_; = ||z;— z;-,|| for ls j <n, 
and requiring second-order continuity with respect to the parameter as shown in 
Fig. 3(b). This chordal parameterization improves on the uniform parameteriza- 
tion of Fig. 3(a), but the splines that we shall develop still have more gentle 
curvature in this case as shown in Fig. 3(c). 

Figure 3 points out the difference between geometric and parametric continu- 
ity. Natural cubic splines have second-order parametric continuity because they 
are described by a function z(t) that is C? continuous. Geometric continuity is a 
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relaxation of parametric continuity: A spline function z(t) has k th-order geomet- 
ric continuity if there is a continuous, monotonic increasing function f such that 
z(f(t)) is C* continuous. In other words, geometric continuity is independent of 
the parameterization. 

Requiring first- and second-order continuity with respect to the parameter 
uses up four degrees of freedom per knot, enough to completely determine a 
parametric cubic spline as shown in Fig. 3(a). In Fig. 3(b), one of these degrees of 
freedom is reclaimed and put to better use by altering the parameter spacing, but 
another degree of freedom can be made available by requiring only continuity of 
slope and curvature. 

Other results on geometric continuity in cubic splines can be found in [1] and 
[2] where Barsky and Beatty show how two extra degrees of freedom can be 
obtained for B-splines by requiring only geometric continuity. We need to obtain 
similar degrees of freedom for interpolating splines, but we shall first concentrate 
on finding good defaults to work from. The new parameters will be of little value 
if it is difficult to set them so as to obtain reasonable results. We cannot afford to 
assume arbitrarily that the best defaults are those that yield some form of 
parametric C? continuity. 

In [8], Manning takes an interesting approach to this problem. He defines a 
specific family of curves so that there is a unique one for each pair of initial and 
final points and directions. He then selects spline directions at each knot so as to 
achieve geometric continuity. His approach provides a certain degree of locality in 
that effects of local perturbations do not propagate past knots where a direction is 
given. 

With Manning’s approach, both degrees of freedom are available to control 
the shape of the curve, and defaults can be selected so as to obtain the most 
pleasing curves. Section 2 explains how to select the defaults by choosing two 
functions and using them to determine the magnitudes of the velocities at each 
knot in such a way as to guarantee that the curves generated will be independent 
of scaling, rotation, and reflection. We can then provide two “tension” parame- 
ters for each knot by simply dividing them into these functions. Essentially the 
same approach would work for other kinds of curves, although there may be more 
parameters to choose. We select parametric cubics here because they are essen- 
tially the simplest curves that can pass through two arbitrary points in two 
arbitrary directions. 

On apparent disadvantage to this approach is the difficulty in solving for the 
directions that provide continuity of curvature. Manning proposes an iterative 
approximation scheme that seems to work well in practice, but he admits that 
there is not always a unique solution and there is no guarantee that the iteration 
always converges to the desired solution. Cubic splines often have very low 
curvature at their endpoints when they have very sharp bends internally, and this 
can introduce extraneous solutions as shown in Fig. 4. The three curves shown are 
all curvature continuous open curves that have given directions at z) and z,, but 
regardless of the initial conditions, Manning’s iteration always converges to one 
of the asymmetrical ones with sharp bends. If zọ is raised and z, lowered until 
the angle z)z,z, is about 122°, the asymmetrical solutions merge with the 
symmetrical ones and the rate of convergence for Manning’s iteration approaches 
zero. 
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z2 


Fig. 4. Three splines of the type proposed in {8}. 


While these kinds of problems do not seem to occur when the angles involved 
are not so large, much additional testing would be necessary in order to verify 
this. In Section 3, we show how all such problems can be avoided by setting up a 
system of linear equations that are easy to solve and guarantee approximate 
continuity of curvature. We derive the specific equations appropriate for the 
families of curves discussed in Section 2, but similar equations could be derived 
for many different classes of curves. 


2. Choosing a Family of Parametric Cubics 


The subproblem to be solved in this section can be stated as follows: Given two 
points z; = (x;, y;) and Z;41 = (X;} 1 Yi+1) and two unit vectors w; and w,,,, find 
an aesthetically pleasing parametric cubic z(t) so that z(0)=z,, z(1)=Z;ņ1 
z’(0) = aw,, and z’(1) = Bw,,,, where a and £ are positive real numbers and z’(t) 
is the componentwise derivative (x’(t), y’(t)) of z(t) = (x(t), y(t)). We also wish 
to introduce two “tension” parameters 7, and 7,,, such that pleasing curves will 
be obtained when +,=7,,,=1, and as the tensions approach oo, the curves will 
approach the line segment joining z; to Z;41- 

In order to guarantee that the results are independent of translation, rotation, 
and scaling, we shall begin by finding a function 2(t) such that 


2(0) = (0,0), 21) = (1,0),  2'(0) = =p(8,4)(c0s8,sind), and 

A 1 : 

2’(1) = Fa 7 ONlcos $, —sing), (2) 
where p and ø are positive real functions to be determined later, 8 = arg w; — 
arg(z,.,—Z,), and ọ = arg(z,,,— 2;)—afgw,.,. We refer to p and o as velocity 
parameter functions because they determine the magnitude of the velocity at t = 0 


and at ¢ =1. (Here arg(x, y) is the angle w such that (x, y) is a positive multiple 
of (cos w,sinw).) We then set 


x; XiT X Vi Yi à T 
= (ea eee ee) ` 6) 


It is not hard to see that the parametric cubic satisfying (2) has Bézier control 
points (0,0), (p/37,Xcos 4, sin@), (1 —(0/37;,ı)cos $, (0/37,,ı)sinġ), and (1,0), 
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so that 


2(1) = £0(1-1)°(cos8, sin) +201- 1)(3- sing eing), (1,0). (4) 


It only remains to choose positive functions p(@,¢) and o(,¢) so that p(6, ) = 


o(, 9) = p(— 4, — 9). 
In [8] Manning chooses 


2 
aCe) = T= east peep 
2 
0(6,%) = rrene rO egosi i 


and then empirically selects c = } to obtain the most pleasing family of curves. 
Here we shall attempt to do a systematic analysis of the vast range of possible 
functions to determine whether slightly more complicated velocity parameter 
functions will yield better results. These functions will have to be evaluated only 
once for each segment of the spline curve, and they have a strong influence on the 
final shape. 


2.1. Mathematical Measures of Smoothness 


One common way of evaluating the smoothness of curves is to integrate the 
square of curvature with respect to arc length. For 2 = (£, ĵ) 


(93 


pee e iy 6 
fras is gr + 92? (6) 


This can be simplified somewhat but it still proves to be analytically intractable. 
This is not surprising considering the complex behavior of numerical solutions. 

Equation (6) is exactly the energy function that the curve of least energy 
minimizes, but if we restrict 2 to be the cubic spline (4), we can investigate the 
velocity parameter functions that minimize (6). Actually we should consider the 
smallest local minimum since (6) approaches 0 as p and o approach oo for fixed 8 
and 9. 

Unfortunately, numerical integration of k? ds proves to be slow and impre- 
cise, and it would have to be repeated a large number of times in order to get a 
good idea what the velocity parameter functions p(6,@) and o(@,¢) should be. 
Instead we shall introduce two other measures of smoothness that behave 
similarly: 


max, |&(1)| (7 


al: (8) 


O<r<l ds 
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Fig. 5. The functions p(@,) and 0(9, @) versus 8 for = 90° — 6, minimizing (a) the magnitude of 
curvature and (b) the magnitude of curvature change. 


Using (7) to measure smoothness corresponds to taking the oo-norm of curvature 
instead of the 2-norm; using (8) gives roughly similar results but applies a greater 
penalty to short periods of relatively high curvature. The velocity parameter 
functions that minimize (8) turn out to be somewhat better behaved than those 
that minimize (7), but the overall character is the same for both measures of 
smoothness. 

For fixed 0 and ¢, the smoothness measures can have multiple local minima 
and the relative smoothness at the local minima can change as 8 and ġo change. 
Therefore, it should not be surprising that the “optimum” velocity parameter 
functions have large discontinuities where they catastrophically change from one 
local minimum to another. When this happens there tend to be (p,o) values 
between the two minima that also generate relatively “smooth” curves, so it is not 
really necessary to use discontinuous functions for p and a. 

Figure 5 illustrates the most basic catastrophe. Near @ = ¢ = 45°, the “opti- 
mum” p increases and the o decreases as 6 decreases. This action tends to reduce 
the curvature where it is maximum near ¢ =1 without introducing other points of 
high curvature. When 0 = @ <4, the situation is entirely different. The high 
curvature near t=1 is best controlled by making o large, and increasing p 
beyond what is needed to control the curvature at t = 0 just makes the problem 
worse. 

When p and o are chosen to minimize (7) as shown in Fig. S(a), there are 
actually two catastrophes, and the short segment between them is particularly 
interesting. Ordinarily, extremely small g values lead to high curvature near ¢ =1, 
but at @=34.1°, ¢=55.9°, the curvature is actually minimized by choosing 
o = 0.261. The choice of p = 2.423 here is extremely critical. As shown in Fig. 6, 
this has the effect of making the last three Bézier control points almost collinear, 
so that the endpoint curvature is not too large in spite of the low velocity. (The 
bold dots in the figure are the Bézier control points of (4).) 


a Fig. 6. The Bézier control polygon for (4) with 9 = 34.1°, @ = 55.9°, where p and 
g are chosen to minimize curvature. 
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3.0 T 


(a) 


Fig. 7. (a) “Optimum” p and ø versus @ when (a) 8 =¢ and (b) 9 = —ẹ for both smoothness 
measures. 


The bizarre situation shown in Fig. 6 does not occur when minimizing the 
derivative of curvature, but there is still a catastrophe near 25°. When 8 + ¢ = 90° 
and @<25° the “optimum” cubic has a point of inflection, but when @ > 25° it 
has none. 

Figure 7 shows how the “optimal” velocity parameter functions grow as 6 
and @ increase. Minimizing (8) produces a catastrophe at 68° where p and o 
increase to about 1.58 and the cubic acquires a single point of maximum 
curvature at ¢ = 0.5. For p = o ~1.34, the maximum curvature occurs at t = 0.08 
and z = 0.92. Intermediate values for p and o produce relatively high change in 
curvature near tf = 0 and ¢ =1. 

Minimizing (7) instead of (8) avoids the catastrophe at 0 = = 68°, but then 
p and o do not approach unique limits as (8, p) > (0,0). Along 0 = — ¢ the limit 
is ¥6 /2, while p and o approach 1 as 8 > 0 when 6 = 9. Under the approxima- 
tion k = d*y/dx*, when either (6) or (8) is used as the measure of smoothness, it 
can be shown that the optimum curves are cubics where pcos@ = ocos¢=1. 
Thus, it seems reasonable to let (p, o) —> (1,1) as (8, ¢) — (0,0). 


2.2, Velocity Parameter Functions 


We are now ready to choose velocity parameter functions p(8,$) and 0(8,¢) so 
as to obtain an aesthetically pleasing family of curves from (4). In order to be 
useful as splines these curves should be fairly easy to compute, and they should 
respond predictably to changes in 6 and ¢. The “optimum” velocity parameter 
functions illustrated in Figs. 5 and 7 fail on both counts, but they still provide 
useful guidelines. The catastrophes in the “optimum” velocity parameter func- 
tions are not features to be preserved, but rather points where many different 
possible values for p and o are feasible. For instance when 0 = 28° and ¢ = 62°, 
good results are obtained for almost any values of p and o as long as p + o = 2.5. 

The actual choice of velocity parameters is necessarily somewhat arbitrary, 
and there is a trade-off between “smoothness” and simplicity. Other properties 
such as approximate extensibility and predictable response to changes are also 
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important, but empirical studies indicate that these goals also tend to be served 
by maximizing “smoothness” and smoothing out the catastrophes. 

We have already decided that p(0,0) = o(0,0) =1. Thus for small angles we 
can approximate the behavior of the curve of least energy and achieve very good 
approximate extensibility. For 0 =, we should approximate the behavior of the 
functions shown in Fig. 7(a), but these functions increase too rapidly or large 
angles: They seem to approach oo well before @ and @ reach 180°. It is 
convenient to let 


2 
= T+¥cos6 forð = > (9) 


p=o 
so as to obtain good approximations to circles. 

Because of the symmetry requirements p(8, $) = 0(¢, 8) = p(— 9, — 9), it 
suffices to choose p and ø for 0 < |8| < $< 180°. Figure 8(a) shows the p(@, >) 
and 0(8, ) that minimize (8) for ¢ = 80°. Figure 8(b) shows practical functions p 
and o that smooth out the catastrophes and are consistent with (9). Similar plots 
for smaller @ would have p and o closer to each other and closer to 1. (The slope 
discontinuities at 64.7° and 69.6° are due to changes in the relative sizes of 
extrema in dk /ds at different parts of the cubic curve.) 

If complexity is of no concern, we might want to choose p(6,@) and o(8, >) 
as follows for 0 < |@| << with angles measured in radians: 


p = f(6,¢)+ (6)sin{ ¥45 


1i 


o 


100,4) - 1(6)sin{ y5), 


aa? +a+2c 
F04) = ar ecasB te’ 


a= -o| a= SS 
y(¢) = ie — 0.15sin(29), 
y(x) = a(x+(x2—1((0.32- 3) mE ral (10) 


A least-squares fit of f to }(p +ø) with p and o chosen to minimize (8) yielded 
a = 0.2678306, c= 0.2638750, d=1.402539, and e = 0.7539063. A possible re- 
finement is to require p < 1.5sin¢/sin(@ + ¢) when ¢ > 47 and —¢<@<0Osoas 
to avoid any possibility of generating a curve with a cusp in it. (This only affects 
the above functions when > 145°.) 

It is desirable to have a simpler approximation that does not use transcenden- 
tal functions other than sines and cosines of @ and . One such approximation is 
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Fig. 8. (a) “Optimum” and (b) practical functions p(8,@) and o(8,¢) versus @ for ọ = 80°. 


the following functions which were developed for the new METAFONT system [6]: 


S ied ee 
P = T4(1— c)cos6 + coos’ 
leaner: a. See 
° = 1+(1=c)cos@ + ccos8’ 
where 
a = a(sinf — bsing)(sing — b sin )(cos 8 —cos@). (11) 


The constants a, b, and c were chosen to minimize an error function based on the 
value of (8) for curves generated by (11) for 116 different (8, p) pairs. This 
suggested a =1.597, b= 0.0700, and c=0.370, but since empirical evidence 
indicated that large values for |p — o| were causing problems, METAFONT uses the 
slightly perturbed values a = V2, b= 3, and c= 4(3- 75). 

Figures 9(a) and (b) show how velocity parameter functions compare in terms 
of the maximum magnitude of change in curvature. The vertical axis in each 
graph gives the ratio of the value of (8) for the indicated velocity parameter 
functions to that for the “optimum” velocity parameter functions. In other words 
for each value of @ and ¢, the graphs give the ratio of (8) to its minimum possible 
value. Since (10) and (11) were derived in order to minimize (8), it is not 
surprising that they perform better than (5) in Figs. 9(a) and (b). We would also 
expect the more complicated functions defined in (10) to perform better than the 
Simpler versions of (11), and this is indeed the case except near 0 = —20°, 
$ = 80°, where /(6, ) is a little too small. 

Figures 10(a) and (b) correspond to Figs. 9(a) and (b), except that the ratios 
are based on the maximum magnitude of curvature instead of the maximum 
magnitude of change in curvature. The graphs show that the curves generated by 
(10) and (11) do have less curvature than Manning’s curves, especially when 9 is 
fairly large and 6 <0. The figures also show that (7) and (8) do behave similarly 
in that they both produce the same overall ranking: Curves from (10) are 
smoother than those from (11) which are smoother than Manning’s curves. 
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Fig. 9. (a), (b) Smoothness ratios based on (8), comparing optimum velocity parameter functions to 
those of (5), (10), and (11). 


Fig. 10. (a), (b) Smoothness ratios based on (7), comparing optimum velocity parameter functions to 
those of (5), (10), and (11). 


Figure 11 shows some of the curves generated by (10) and (11). They are 
similar for moderate angles, but the simpler equations set p too small and o too 
large when @ = — 90°. Equation (11) does not perform well in such extreme cases 
because it does not allow p—o to be large enough when 0 < @/¢ <1 without 
making o — p too large when —1 < 8/¢ <0 or moving the crossover point where 
p =o too close to 6/¢ = 0. 


(a) (b) 


Fig. 11. (a) Curves from (10) and (b) curves from (11), both with @ = 45°. 
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3. Mock Curvature Constraints 


Here we need to extend the notation of Section 2 so that 8, = arg w, — arg(z,,, — Z,) 
for O< j<n, $,=arg(z,—z,_,)—argw, for 0< j<n, and d, is the Euclidean 
length of the vector z,,,;—z; for O< j<n. If the problem is to find a closed 
curve with no directions given, it will be convenient to sometimes use alternative 
names Z,.1, 9,, O15 Tp and 7,,, for z,, Oos $1, To, and 7,, respectively. We can 
then define y, = arg(z;,,— z,)—arg(z;— 2,1) for 1< j <n’, where n’=n for 
closed curve problems with no directions given and n’ =n —1 otherwise. Unless 
stated otherwise, all y, are at most 180° in absolute value. 

Where w, has been given in advance, it simply determines ġ; and 6; other w, 
need to be determined by solving for 6, and ¢,. Since the problem of finding 
direction angles can be broken into independent subproblems separated by knots 
where directions are given, we can assume that no directions other than wọ and w, 
are given. For closed curve problems we can assume that no directions at all are 
given, otherwise the problem could be reduced to one or more open curve 
problems. 

The requirement that the curvature be continuous at some knot z,,0<i<n, 
is 

ky (2,215 25% 1s GT) = kozi Zap Wo Wan To T41) 
where kọ and k, are functions that give the curvature at ¢=0 and ¢=1 in terms 
of the endpoints, terminal directions, and tension parameters for the family of 


curves being used. Because of the requirement for invariance under translation, 
rotation, and scaling, there exists a function k such that 


kl Zj Zj Wo Mar To Tai) = kh8, b1 T T41)/4, 
and 
kalZ, Zj Mar TT) = kloan b F415)/4,- (12) 


Any particular family of curves determines a specific function k that satisfies 
(12). The corresponding mock curvature function k consists of the linear terms in 
the Taylor series for k(0,$,7, 7), expanded about (0, ¢) = (0,0). For the curves 
determined by (3) and (4) with p and o determined by (10) or (11), 
20(8,)sin(6 + o)/7 —6sin8 


k(0,6,7,7) = 5 
TE (o(0,4)/7) 


and 
Ue aoe a (2229) 66), (13) 


where the angles are measured in radians. Since the tension parameters are always 
known in advance, they are treated like constants in this expansion. 
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Instead of requiring continuity of curvature, we will obtain a system of linear 
equations by requiring the mock curvatures to match at each knot. When the 
angles 6, and 9, are small, this guarantees that the resulting splines will have 
approximate second-order continuity. For instance 


|k(8,,1,1)— k(8,4,1,1)| 


eo <. 

max(|k(6,¢,1,1) |, |K(, 4,1, 1) |) a 

when |8], |p| < 22.5° and p and o are determined by (11) with a= V2, b= k, 
and c= }(3— V5). 

In many applications the angles 0; and ¢, seldom get much larger than 22.5° 
but experience has shown that the splines usually appear very smooth even when 
some @, and 9, are much larger than 22.5°. An iterative algorithm such as 
Manning’s can be used to obtain true curvature continuity if desired, but this 
necessitates giving up the simplicity and guaranteed behavior of a system of linear 
equations. 

Continuity of mock curvature requires 


Ko Oavirtiaa) Hitomi) 9 wercicw. (4 


di- 
Combining this with the first-order continuity equations 
64+6,=-y, frlsi<n (15) 


gives enough equations to determine all 0, and ¢, for closed curve problems. For 
open curve problems when directions wọ and w, are given in advance, these 
provide the necessary additional equations by fixing 4 and ¢,, but otherwise 
additional constraints are needed. 

The additional constraints are controlled by special “curl” parameters x and 
Xa- There should be one such constraint for each endpoint where no direction is 
specified. They have the form 


k (8,1, 771) a Xok ($1, bos F To) (16a) 
and 
K( bys p15 Tas M1) ad X nk (0,15 bn Ta- Ta) (16b) 


The curl parameters give the ratio of the mock curvature at the endpoints to that 
at the adjacent knots. They should probably have default values of 1 so that the 
first and last spline segments will usually be good approximations to circular arcs 
as in Fig. 12(b). 

We now have a system of equations consisting of (14), (15), and possibly 
(16a) and/or (16b). If 6 or ¢, have been given in advance then they may be 
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zo 22 Zo 22 20 z 
(a) (b) (c) 
Fig. 12. (a) Xo = Xz = 0. (©) Xo = X2 =1. © Xo = X2 = ©. 


regarded as constants. The first step is to rewrite (16a) and (16b) as 


3 + xoF?(3%) —1) T + x,7°.,(37,-1) 
8. = 0 oTi 0 d = n n'n-1 n 6 17 
s To 3R ~1)+ x07? os n P Brna Tl) Xni a ( 


so that 8 and ¢, can be eliminated. Then (15) may be used to eliminate all ¢, so 
that the remaining variables are 6,,6,,...,6,,, and the remaining equations are 
those given by (14) with appropriate substitutions. This system has some im- 
portant properties that may be summarized as follows. 


Theorem 1. If n> 2, if all tension parameters satisfy the bound 1,,7,> Tin > 3, 
and if any curl parameters satisfy Xo, X , = 0, then after the aforementioned substitu- 
tions, all coefficients of 8,,0,,...,6,, in (14) are nonnegative, and for each i, the 
coefficient of 0, is at least 3T pin ~ 1 times the sum of all the other coefficients in that 


equation. 


Proof. The bounds on the tension parameters guarantee that the coefficient of 8 
in (13) will be negative, the coefficient of ọ will be positive, and the magnitude of 
the former will be at least 37 „in — 1 times the latter. When 1<i<n’ in (14), the 
my relevant substitutions are $, =—y,— ġ; for |j — i| <1, so the coefficients of 

1 9,, and 6,,, clearly have the required properties. For closed curve problems, 
fe same holds for i=1 and i = n’, otherwise additional substitutions eliminate b 
and @¢, so that klon bos Tis To) depends only on ġġ, and k(O, -i Ons T—19 Tr) 
depends only on ĝ„_;ı. We need only show that both of these variables have 
nonpositive coefficients. This is clearly true for given directions, and it also holds 
when curl parameters apply since the coefficients in (17) are at most 37) —1 and 
37, — 1, respectively. o 

Theorem 1 shows that subject to certain reasonable limitations on the tension 
and curl parameters, the system of equations is diagonally dominant, and hence it 
has a unique solution. Actually the solution is unique only up to the choice of the 
angles y,. Ordinarily all ~, should be chosen so that they are at most 180° in 
absolute value, but it is possible to add multiples of 360° to them. The effect of 
such a change is usually to add a loop to the curve as in Fig. 13. 

Theorem 1 also shows that the splines have approximate locality in the sense 
that changes in direction angles fall off exponentially as one moves away from a 
disturbance. Specifically, suppose a given direction 9, is displaced by an angle ô 
and let A be the matrix of coefficients of 0,,0,,...0,, from (14) after the 
substitutions. The change in 6,, 8,...,6,, due to this displacement is given by the 
solution vector © to A@ = Se, where e; = (1,0,0,...,0)”. 
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22 
22 


20 z 3 24 
20 21 23 24 


(a) (b) 
Fig. 13. A spline computed (a) with y, = —90° and (b) with y, = 270°. 


K .  Fig.14. Exponential decline in the effect 
zo #1 22 23 74 of a 45° change in direction. 


We know that A is tridiagonal with nonnegative entries, and within each row 
the diagonal element dominates the sum of the other two elements by at least a 
factor of 37,,,, — 1. It is not hard to see that for any two adjacent components of 
O, either ©, = 0 or 8, _,/O, $1—3t,;,. This is trivial for k =n’, and it may be 
extended inductively to smaller k using the fact that A,, 2 (3Tmin ~ DC Ag. 4-1 + 
Ax 441): Thus, J knots away from where a given direction is changed; the effect 
of the change is reduced by at least a factor of (3Tmin— 1)”. In practice the 
reduction is often by a somewhat greater amount as in Fig. 14 where Tmin = 1 and 
9/8, = — fi- 

When a knot z, is displaced, three mock curvature constraints are directly 
affected due to changes i in d,- di Y,-1, Yi and p,,,. The adjustment will cause 
some change in ¢,_, and 6,,,, and the effect on 6,,, and 6,_, for j>l is 
equivalent to what would happen if directions w,,, and w,_, were given in 
advance. The change in 6,,, will be at most 1/(3t;,—1)/~* as great as the 
change in 8,,,, and the change i in 6,_, will be at most 1/(3fmin — 1)/~* as great as 
the change in ¢,_,. If the original ‘problem were to find a ‘lobed curve with no 
directions given, then these two effects will add together so that the change in 6, ,, 
will be at most 1/(37,,;, ~1)/~! of the change in @,,, plus 1/(3Tmin —1)"~'/ of 
the change in ¢,_. 


Tmin 


4. Conclusion 


We have developed a tridiagonal system of linear equations that can be solved in 
linear time to determine the spline direction at each knot so as to match mock 
curvatures. It is necessary to use arctangents to set up the system of equations, 
and to use sines and cosines to recover the resulting spline directions; but this 
work can be reduced to one arctangent, one sine, and one cosine per knot on the 
spline. When all unit direction vectors w, have been determined, the ith spline 
segment is the cubic whose Bézier control points are 


q Peisa) gte ofhi) 


Zi Zi [z1 — 2M, Zisi 371 2:41 — zilwis 
i 


and Z;4, 
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(a) (b) 
ö 3 4 ð 3 4 
1 2 2 1 2 
(d) (e) (fy 


Fig. 15. (a) 7,=%=0.75. (b) 1=%=0.85. (c) 7 =F=1. (d) 7 =F =12. (e) =F =15. (f) 


naa. 


Thus, the total work required in order to find a computationally convenient 
representation of the desired spline is linear in the number of knots. The constant 
factor obviously depends on the implementation, but it is usually not much higher 
than the similar factor for natural cubic splines. 

The splines do not have locality, but we have shown that changes in direction 
angles fall off exponentially. The rate of decline depends on how small the 
tension parameters are allowed to be, but at least a factor of 2 per knot is 
guaranteed for the default tensions, and decay rates as high as a factor of 4 per 
knot are typical. It should be noted that an exponential decline in angular change 
does not guarantee that curve displacements decline similarly because it is 
technically feasible for d, to be exponential in j. 

The curve families discussed in Section 2 and defined by (10) and (11) are 
somewhat arbitrary, and the concept of mock curvature could be applied to other 
families of curves. It would be desirable to find p and o functions simpler than 
(10) that perform better than (11), although even the simplified functions of (11) 
produce very good results for problems such as that shown in Figs. 2(a) and 3(c). 
In fact the splines produced by (11) have performed outstandingly in extensive 
tests with Knuth’s new METAFONT system [6]. It is easy to compute the splines by 
solving the equations given in Section 3, but readers wishing a detailed descrip- 
tion of one implementation can refer to [6]. 

The tension parameters have proved very useful in METAFONT, because it is 
often much easier to control a shape by specifying tensions than by specifying 
directions or giving more knots. Figure 15 shows the effect of adjusting all the 
tension parameters simultaneously, and Fig. 16 shows the effect of adjusting only 
one tension parameter and leaving the rest of the tension parameters at the 
default values (7, = 7,=1). Note how the tension effects the direction chosen at 
knot 2 in Fig. 16. 
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(a) 
Fig. 16. (a) 3 =1.2. (b) % =1.5. (0) % =2. 


Acknowledgments 


The author is indebted to Professor Donald Knuth of Stanford University for doing the first complete 
implementation and helping to refine some of the ideas. Professor Kevin Karplus of Cornell 
University helped with the numerical investigation of the “optimal” p and o functions. 


References 


1. Brian A. Barsky, The Beta Spline: A Local Representation Based on Shape Parameters and 
Fundamental Geometric Measures, Ph.D. thesis, Univ. of Utah, December, 1981. 

2. Brian A. Barsky, John C. Beatty, Varying the betas in beta-splines, Univ. of California, Berkeley 
TR CSD 82/112, December 1982. 

3. Tony D. DeRose and Brian A. Barsky, Geometric continuity and shape parameters for Catmull- 
Rom Splines, Graphics Interface ’84, 57-64. 

4. M. P. Epstein, Parametric interpolation, SIAM Journal of Numerical Analysis 13 (1976), 261-268. 

5. B. K. P. Horn, The curve of least energy, Massachusetts Institute of Technology A.I. Memo No. 
612, January 1981. 

6. Donald E. Knuth, Computers and Typesetting, Vol. 4: METAFONT the Program, Addison 
Wesley, to appear. 

7. Doris H. U. Kochanek and Richard H. Bartels, Interpolating splines with local tension, continu- 
ity, and bias control, Computer Graphics 18 (1984), 33-41. 

8. J. R. Manning, Continuity conditions for spline curves, Computer Journal 17 (1974), 181~186. 

9. Even Mehlum, Curve and Surface Fitting Based on Variational Criterie for Smoothness, Oslo, 
1969. 


Received January 31, 1985, and in revised form September 4, 1985. 


